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Abstract 


We  have  explored  the  degree  to  which  an  intermolecular  potential  for  the  explosive 
hexahydro- 1 ,3,5-trinitro-l  ,3,5-s-triazine  (RDX)  (J.  Phys.  Chem.  vol  101B,  p.  798)  is  transferable 
for  predictions  of  crystal  structures  (within  the  approximation  of  rigid  molecules)  of  a  similar 
chemical  system,  in  this  case,  polymorphic  phases  of  the  2,4,6,8,10,12- 
Hexanitrohexaazaisowurtzitane  (HNIW)  crystal.  Molecular  packing  and  isothermal-isobaric 
molecular  dynamics  calculations  performed  with  this  potential  reproduce  the  main 
crystallographic  feature  of  the  £-,  P-and  y-HNIW  crystals.  Thermal  expansion  coefficients 
calculated  using  the  present  model  predict  near  isotropic  expansion  for  the  e-  and  Y-HNIW 
crystals  phases  and  anisotropic  expansion  for  P-HNIW. 
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1.  Introduction 


Satisfying  ever-increasing  demands  for  inexpensive,  efficient,  and  rapid  development  of 
high-performing  energetic  materials  is  extremely  challenging.  Synthesis  and  testing  of  energetics 
is  often  costly,  dangerous,  and  time  consuming;  thus,  a  screening  mechanism  is  needed  to  reduce 
unnecessary  measurements  on  unsuitable  candidate  materials.  Such  screening  can  be  accomplished 
through  the  application  of  theoretical  chemical  models.  These  models  can  be  used  to  predict  simple 
thermodynamic  properties  of  die  candidates,  such  as  heats  of  formation  or  the  density  of  the  material, 
which  are  sufficient  to  aid  formulators  in  decisions  of  whether  a  candidate  energetic  material 
warrants  further  investigation.  Due  to  timely  advances  in  parallel  computer  architectures,  there  is 
great  promise  that  atomistic  modeling  will  become  integral  to  the  development  process.  In  addition 
to  being  efficient  screening  tools,  theoretical  chemical  models  can  also  provide  atomic-level 
information  that  could  not  otherwise  be  readily  obtained  through  measurement. 

Our  efforts  at  developing  accurate  predictive  models  of  energetic  crystals  have  included  the 
development  of  an  intermolecular  potential  that  describes  the  structure  of  the  a-form  of  the 
hexahydro-1 ,3,5,-trinitro-l  ,3,5-s-triazine  (RDX)  crystal,  one  of  the  most  commonly  used 
explosives  [1].  The  potential  energy  function  that  describes  the  system  is  composed  of  pairwise 
atom-atom  (6-exp)  Buckingham  interactions  with  explicit  inclusion  of  the  electrostatic  interactions 
between  the  charges  associated  with  various  atoms  of  different  molecules.  The  parameterization  of 
the  potential  function  was  done  such  that  molecular  packing  (MP)  calculations  reproduce  the 
experimental  structure  of  the  crystal  and  its  lattice  energy.  Isothermal-isobaric  molecular  dynamics 
(NPT-MD)  simulations  using  this  potential  energy  function  predicted  crystal  structures  in  excellent 
agreement  with  experiment.  The  main  limitation  of  the  model  is  due  to  the  assumption  of  rigid 
molecules.  Nevertheless,  it  can  be  used  to  study  processes  at  temperatures  and  pressures  where 
molecular  deformations  are  negligible. 

The  development  of  a  simple  model  that  accurately  represents  a  specific  chemical  system  is 
significant.  However,  the  utility  of  the  model  is  substantially  enhanced  if  it  reasonably  represents 
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a  series  of  chemical  systems  rather  than  a  single  one.  This  investigation  explores  the  degree  of 
transferability  of  the  RDX  crystal  potential  energy  function  [1]  to  another  cyclic  nitramine 
crystal — that  is,  the  new  explosive,  2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  (HNIW)  (see 
Figure  1). 

HNIW,  a  polycyclic  nitramine,  has  been  characterized  as  “the  densest  and  most  energetic 
explosive  known”  [2].  It  exists  in  at  least  five  polymorphic  states  [3],  four  of  which  (a-hydrate,  e, 
P,  and  y)  are  stable  at  ambient  conditions  and  have  been  resolved  by  x-ray  diffraction  [4].  The 
a-hydrate  phase  has  an  orthorhombic  structure  with  Pbca  symmetry  and  with  Z  =  8  molecules  per 
unit  cell;*  the  e  polymorph  crystallizes  in  the  P2/n  space  group  and  has  Z  =  4  molecules  per  unit 
cell;  P-HNIW  has  Pb2ja  symmetry,  with  Z  =  4  molecules  per  unit  cell;  and  the  y  phase  has  P2j/n 
symmetry,  with  Z  =  4  molecules  per  unit  cell.  Figure  1  illustrates  these  polymorphs  of  HNIW;  the 
C  polymorph  is  evident  at  high  pressure,  but  its  crystal  structure  has  not  yet  been  resolved  [4].  As 
evident  in  the  figure,  the  molecular  structure  of  the  polymorphs  differ  mainly  in  the  orientation  of 
the  nitro  groups  relative  to  the  ring.  Two  different  rankings  of  the  relative  stabilities  of  these 
polymorphs  have  been  proposed:  a-hydrate  >  e  >  a-anhydrous  >  P  >  y  [5]  and  e  >  y  >  a-hydrate 

>  P  [6]. 


2.  Intermolecular  Interaction  Potential 


For  this  study,  we  use  the  same  values  of  the  6-exp  potential  parameters  as  in  the  RDX  study  [1]. 
The  Coulombic  terms  are  determined  by  fitting  partial  charges  centered  on  each  atom  of  the  HNIW 
molecule  to  a  quantum  mechanically  derived  electrostatic  potential,  performed  at  the  HF/6-31G** 
level  as  described  in  the  earlier  study  [1].  For  the  treatment  of  different  phases  of  the  HNIW  crystal, 
we  also  assume  that  there  exists  a  transferability  of  the  electrostatic  charges.  Consequently,  we  have 
used  in  all  calculations  a  single  set  of  atom-centered  charges  determined  for  an  HNIW  molecule  with 


*  The  x-ray  diffraction  data  correspond  to  die  hydrated  crystal,  in  which  the  number  and  orientation  of  the  water 
molecules  are  not  resolved  [7]. 
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Figure  1.  Molecular  Configuration  of  HNIW.  Atom  Labels  Are  Consistent  With  the  Indices  Given  in  Table  1.  Structures  in  the 
Lower  Portion  of  the  Figure  Illustrate  the  Differences  Between  the  Polymorphic  Phases;  the  Hydrogen  Atoms  Are  Not 
Shown  for  Clarity. 


atomic  arrangement  consistent  with  the  crystallographic  configuration  of  the  a-  polymorph  [4].*  The 
resulting  electrostatic  charges  are  given  in  Table  1,  and  the  parameters  for  the  6-exp  terms  are  given 
in  Sorescu,  Rice,  and  Thompson  [1]. 

2.1  Details  of  the  Calculations.  The  tests  of  this  potential  were  done  by  MP  calculations  with 
and  without  symmetry  constraints  and  NPT-MD  simulations  at  zero  pressure  over  the  temperature 
range  4.2-425  K  for  the  e,  P,  and  y  polymorphs  of  HNIW.  We  do  not  report  results  of  a-HNIW, 
since  no  x-ray  diffraction  data  are  available  for  the  anhydrous  form  of  this  polymorph.  Also,  the 
x-ray  diffraction  data  for  a-hydrate  HNIW  is  incomplete.**  Thus,  a  suitable  comparison  between 
theoretical  predictions  and  experiment  for  a-HNIW  is  not  possible  at  this  time.  In  all  calculations 
reported  in  this  study,  the  crystal  is  represented  as  an  ensemble  of  rigid  molecules.  The  independent 
degrees  of  freedom  are  the  six  unit  cell  constants  (a,  b,  c,  a,  P,  y),  the  three  rotations  (0X,  02, 03),  and 
the  three  translations  (tl5  x2,  t3)  for  every  molecule  considered  in  the  simulation.  Details  of  the  MP 
and  NPT-MD  calculations  are  described  in  Sorescu,  Rice,  and  Thompson  [1]  and  are  briefly 
summarized  in  following  text. 

2.1.1  MP  Calculations.  Two  series  of  MP  calculations  (minimizations  of  the  lattice  energy  with 
respect  to  the  structural  degrees  of  freedom  of  the  crystal)  were  performed.  In  the  first  series,  using 
the  program  PCK91  [8],  the  space-group  symmetries  of  the  crystals  were  maintained  throughout  the 
energy  minimization.  In  these  calculations,  all  three  dimensions  of  the  unit  cell  and  the  three 
rotations  and  translations  of  the  molecule  in  the  asymmetric  unit  were  allowed  to  vary,  except  for 
the  case  of  the  P-phase,  where  the  translation  along  the  b  axis  was  frozen  due  to  the  symmetry 
restrictions.  The  cell  angles  a,  p,  and  y  were  frozen  at  90°  for  the  calculations  of  the  P-polymorph. 
For  the  e-  and  y-polymorphs,  the  angle  P  of  the  unit  cell  was  allowed  to  vary,  and  a  and  y  were 


*  MP  (PCK91)  calculations  using  HF/6-31G**  atomic  charges  determined  for  HNIW  with  molecular  arrangements 
consistent  with  the  p-  and  y-  polymorphic  crystallographic  configurations  predict  lattice  parameters  similar  to  those 
that  use  HF/6-3 1G**  atomic  charges  assuming  the  a-hydrate  HNIW  molecular  configuration.  Cell  edge  lengths  differ 
by  no  more  than  0.015  A,  and  cell  angles  differ  by  less  than  0.17°.  The  electrostatic  energies  differ  by  less  than 
5  kJ/mol. 

The  x-ray  diffraction  data  correspond  to  the  hydrated  crystal,  in  which  the  number  and  orientation  of  the  water 
molecules  are  not  resolved  [7]. 
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Table  1.  Electrostatic  Charges  for  the  HNIW  Molecule 


Atoma 

Charge/|e| 

Atom4 

Charge/ |e| 

Cl 

0.256850 

01 

-0.362942 

C2 

0.183418 

02 

-0.369093 

C3 

0.020184 

03 

-0.395097 

C4 

0.089133 

04 

-0.377983 

C5 

0.548294 

05 

-0.357808 

C6 

0.368468 

06 

-0.431362 

N1 

-0.196698 

07 

-0.395640 

N2 

0.680596 

08 

-0.372372 

N3 

-0.357766 

09 

-0.401700 

N4 

0.786509 

010 

-0.416611 

N5 

-0.326781 

Oil 

-0.381535 

N6 

0.754647 

012 

-0.425342 

N<7 

-0.391804 

HI 

0.112444 

N8 

0.753909 

H2 

0.117717 

N9 

-0.577874 

H3 

0.182501 

0.844354 

H4 

0.201482 

L Nu 

-0.234524 

H5 

0.081152 

0.696675 

H6 

0.094597 

a  The  atom  indices  are  consistent  with  the  labels  in  Figure  1. 


frozen  at  90°.  The  positions  and  orientations  of  all  other  molecules  in  the  unit  cell  are  determined 
through  symmetry  operations  relative  to  the  molecule  in  the  asymmetric  unit. 

A  second  series  of  MP  calculations  were  performed  in  which  the  crystal  symmetries  are  not 
constrained;  the  methods  are  described  in  our  previous  paper  [1],  and  the  algorithm  used  is  that  in 
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the  program  LMIN  [9].  hi  these  calculations,  the  cutoff  parameters  P  and  Q  [1]  were  set  equal  to 
17.5  and  18.0,  respectively. 

2.12  Isothermal-Isobaric  Molecular  Dynamics  Calculations.  Simulations  of  the  e-,  (5-,  and 
y-  phases  of  the  HNIW  crystal  at  various  temperatures  in  the  range  4.2-425  K  and  at  zero  pressure 
were  performed  using  the  algorithm  proposed  by  Nose  and  Klein  [10]  as  implemented  in  the 
program  MDCSPC  [11].  Details  of  the  calculations  are  given  in  Sorescu,  Rice,  and  Thompson  [1] 
and  remain  the  same,  with  the  following  exceptions:  The  MD  simulation  cell  consists  of  a  box 
containing  12  (3x2x2),  27  (3x3x3),  and  12  (2x3x2)  crystallographic  unit  cells  for  phases  e,  P,  and 
y,  respectively.  The  lattice  sums  were  calculated  in  all  dimensions  in  these  simulations.  The 
interactions  were  determined  between  the  sites  in  the  simulation  box  and  the  nearest-image  sites 
within  the  cutoff  distance.  The  cutoff  distances  are  =  10.01  A,  1 1.61  A,  and  9.80  A  for  phases 
e,  P  and  y,  respectively.  For  simulations  at  4.2  K,  the  initial  positions  of  the  atoms  are  identical  to 
those  in  the  experimental  structures.  The  equations  of  motion  were  integrated  for  2,000  time  steps 
(1  time  step  =  2x1  O' 15  s)  to  equilibrate  the  system.  After  the  equilibration,  an  additional  5,000  (for 
temperatures  below  300  K)  or  8,000  (for  temperatures  between  300  and  425  K)  steps  were 
integrated,  during  which  average  properties  were  calculated.  In  subsequent  runs  for  successively 
higher  temperatures,  the  initial  atomic  positions  and  velocities  are  identical  to  those  obtained  at  the 
end  of  the  preceding  lower-temperature  simulation.  The  velocities  were  scaled  over  an  equilibration 
period  of  2,000  steps,  in  order  to  achieve  the  desired  external  temperature  and  pressure,  followed  by 
the  5,000-  or  8,000-step  integration  for  calculation  of  averages.  The  cumulative  mass-center  radial 
distribution  functions  (RDF)  and  averages  are  calculated  for  the  mass  centers  and  Euler  angles  of 
the  molecules.  The  RDFs  and  averages  were  obtained  from  values  calculated  at  every  tenth  step 
during  the  trajectory  integrations. 

2.2  Results  and  Discussions. 

2.2.1  MP  Calculations.  The  results  of  the  MP  calculations  with  (denoted  as  LMIN)  and  without 
(denoted  PCK91)  symmetry  constraints  are  given  in  Table  2.  The  relaxation  of  the  symmetry 
conditions  has  a  very  small  effect  on  the  final  geometric  parameters,  and  both  sets  of  MP 
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Table  2.  Lattice  Parameters  for  Polymorphs  of  HNIWa 


Lattice 

Parameter 

Experiment 

[4] 

PCK91 

LMIN 

NPT-MDb  | 

4.2  K 

300  K  | 

|  e-HNIW 

a  (A) 

8.8278 

8.7973  (-0.4) 

8.7948  (-0.4) 

8.7953  (-0.4) 

8.8420  (0.2) 

b(A) 

12.5166 

12.4986  (-0.1) 

12.4999  (-0.1) 

12.5006  (-0.1) 

12.5837  (0.5) 

c(A) 

13.3499 

13.4071  (-0.4) 

13.4055  (0.4) 

13.4066(0.4) 

13.4897  (1.0) 

a(°) 

90.000 

90.000c 

89.998  (0.0) 

90.000  (0.0) 

89.999  (0.0) 

P  n 

106.752 

105.150  (-1.5) 

105.137  (-1.5) 

105.134  (-1.5) 

105.377  (-1.3) 

Y(°) 

90.000 

90.000° 

90.001  (0.0) 

90.000(0.0) 

90.012  (0.0) 

p-HNIW 

a  (A) 

9.6764 

9.5242  (-1.6) 

9.5239  (-1.6) 

9.5272  (-1.5) 

9.6106  (-0.7) 

b(A) 

13.0063 

12.8726  (-1.0) 

12.8728  (-1.0) 

12.7485  (2.0) 

12.9316  (-0.6) 

c(A) 

11.6493 

11.7025  (0.5) 

11.7020(0.4) 

11.7050(0.5) 

11.7566(0.9) 

a(°) 

90.000 

90.000° 

90.000  (0.0) 

89.999  (0.0) 

90.009  (0.0) 

P(°) 

90.000 

90.000° 

89.999  (0.0) 

90.000  (0.0) 

90.009  (0.0) 

Y(°) 

90.000 

90.000° 

90.000  (0.0) 

90.000(0.0) 

89.999  (0.0) 

j  y-hniw  1 

a  (A) 

13.2310 

13.4342  (1.5) 

13.4348  (1.5) 

13.4370  (-1.6) 

13.5144(2.1) 

b  (A) 

8.1700 

7.9095  (-3.2) 

7.9074  (-3.2) 

7.9092  (-3.2) 

7.9690  (-2.5) 

c(A) 

14.8760 

14.8531  (-0.2) 

14.8519  (-0.2) 

14.8548  (-0.1) 

14.9413  (0.4) 

cc(°) 

90.000 

90.000° 

90.000  (0.0) 

90.000(0.0) 

89.986  (0.0) 

P(°) 

109.170 

108.84(0.3) 

108.734  (-0.4) 

108.863  (-0.3) 

109.064  (-0.1) 

Y(°) 

90.000 

90.000° 

90.001  (0.0) 

90.000  (0.0) 

90.013  (0.0) 

a  Percent  deviations  from  experiment  in  parentheses. 
b  Time-averaged  values  (see  text). 
c  Fixed  throughout  calculation. 
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calculations  predict  almost  identical  structures.  Deviations  of  the  cell  lengths  from  experiment  are 

less  than  0.4%,  1.6%,  and  3.2%  for  the  e,  P,  and  y  phases,  respectively.  For  those  crystal  symmetries 

in  which  the  P  cell  angle  was  allowed  to  vary  during  minimization  (the  e  and  y  phases),  deviations 

of  this  angle  from  experiment  are  no  greater  than  1.5%.  Lattice  energies  per  molecule  for  the  e,  P, 

and  y-HNIW  phases  (-210.47  kJ/mol,  -207.43  kJ/mol,  and  -201.42  kJ/mol,  respectively)  support 

the  polymorph  stability  ranking  of  e  >  P  >  y  given  by  Russell  et  al.  [5].  Small  differences  in  the  total 

lattice  energies  (<0.5  kJ/mol)  between  the  constrained  and  unconstrained  calculations  are  due  to  the 

fact  that  the  unconstrained  simulations  do  not  use  the  accelerated  convergence  method  for  evaluation 

of  the  —  lattice  sums. 
r6 

2.2.2  NPT  Molecular  Dynamics  Calculations.  The  analysis  of  time  histories  (not  shown)  of 
the  lattice  parameters  (a,  b,  c,  a,  P,  and  y)  indicate  that  these  properties  are  well  behaved  after  the 
equilibrium  is  reached  (i.e.,  each  parameter  oscillates  about  die  average  value  for  the  duration  of  the 
trajectory).  Time  histories  of  the  rotational  and  translational  temperatures  and  for  the  pressure  show 
similar  behavior,  indicating  that  the  system  is  at  thermodynamic  equilibrium. 

The  crystal  structure  information  resulting  from  NPT-MD  simulations  at  zero  pressure,  T  =  4.2 
and  300  K,  is  given  in  Table  2.  The  lattice  dimensions  obtained  at  T  =  4.2  K  are  in  very  close 
agreement  with  those  determined  in  the  MP  calculations.  This  is  expected,  since  the  thermal  effects 
at  4.2  K  should  be  minimal  and  the  thermal  averages  at  this  temperature  should  be  close  to  the  values 
corresponding  to  the  potential  energy  minimum.  At  300  K,  the  average  lattice  dimensions  agree  very 
well  with  the  experimental  values — the  corresponding  differences  for  a,  b,  and  c  cell  lengths  being, 
respectively,  0.2%,  0.5%,  and  1.0%  for  the  e  phase;  0.7%,  0.6%  and  0.9%  for  the  P  phase;  and  2.1%, 
2.5%,  and  0.4%  for  the  y  phase.  For  the  e  and  y  phases,  the  variations  of  the  unit  cell  angle  P  from 
the  experimental  values  are  1.3%  and  0.1%,  respectively,  while  the  other  two  angles  of  the  unit  cell 
remain  approximately  equal  to  90°.  For  the  P  phase,  all  three  crystallographic  angles  remain 
approximately  equal  to  90®,  in  agreement  with  experiment. 

Figure  2  provides  a  visual  comparison  of  the  average  mass-center  fractionals  and  Euler  angles 
for  each  of  the  four  molecules  within  the  unit  cell  of  e-HNIW  with  experimental  values;  the  data  for 
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Figure  2.  Comparison  of  the  Time-Averaged  Mass  Center  Fractional  Positions  and  Euler 
Angles  (X-Convention  [12])  With  Experimental  Results  at  300  K  and  0  atm  for  the 
e  Phase.  These  Time  Averages  for  Each  Molecule  in  the  Unit  Cell  Are  Over  All  Unit 
Cells  in  the  Simulation  Box. 


9 


the  three  polymorphs  of  HNIW  are  given  in  Tables  3  and  4.  Increasing  the  temperature  from  4.2  to 
300  K  does  not  produce  any  significant  displacement  of  the  molecular  mass-centers  or  increase  the 
degree  of  rotational  disorder.  In  addition,  there  is  a  slightly  better  agreement  with  experiment  for 
the  orientational  parameters  and  fractional  coordinates  of  the  molecules  at  300  K  than  at  lower 
temperatures. 

Similar  good  agreement  with  experiment  exists  for  predictions  of  0-  and  y-HNIW  (not  shown). 
The  largest  deviation  between  experiment  and  predictions  of  molecular  orientation  occurs  for  the 
Euler  angle  $  for  y-HNIW;  the  predicted  value  is  4.4°  larger  than  the  experimental  value. 

The  mass-center-mass-center  RDFs  for  the  c  polymorph  (Figure  3)  exhibit  well-ordered 
structures  with  correlation  at  long  distances  at  higher  temperatures.  The  positions  of  the  major  peaks 
do  not  change  significantly,  and  the  main  temperature  effect  is  the  broadening  of  the  peaks  and 
partial  overlapping  of  some  of  them.  The  features  of  the  RDFs  for  the  P  and  y  phases  are  similar 
to  those  shown  in  Figure  3. 

Finally,  the  cell  edge  lengths  and  volumes  have  a  linear  dependence  on  temperature  over  the 
range  4.2-425  K,  as  evident  in  Figure  4.  Linear  and  volume  thermal  expansion  coefficients  are 
determined  from  linear  least-squares  fits  to  these  data,  the  parameters  of  which  are  listed  in  Table  5. 
The  linear  and  volume  expansion  coefficients, 


(1) 


where  X  denotes  cell  edge  lengths  a,  b,  or  c  or  volume  V  and  Cx  is  the  corresponding  thermal 
expansion  coefficient,  were  calculated  at  T  =  300  K;  they  are  given  in  Table  5  for  the  three 
polymorphic  phases.  The  results  indicate  a  near  isotropic  thermal  expansion  for  the  e  and  y  phases 
(coefficients  are  within  ~20%  of  one  another),  and  anisotropic  expansion  for  p-HNIW.  For 
P-HNIW,  the  linear  expansion  coefficient  for  cell  edge  a  is  ~70%  larger  than  those  for  b  and  c.  The 
volume  expansion  coefficients  for  the  three  polymorphs  have  similar  values  (~6  x  10‘5  K'1).  At 
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Table  3.  Average  Values  of  the  Fractional  Coordinates  and  Euler  Angles  (X-Conventiona)  of  the  Molecular  Mass  Center 
Calculated  From  Trajectories  at  4.2  and  300  K 
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Table  4.  Average  Values  of  the  Fractional  Coordinates  and  Euler  Angles  (X-Conventiona)  of  the  Molecular  Mass  Center 
Calculated  From  Trajectories  at  4.2  and  300  K  With  the  Corresponding  Values  Determined  for  the  Experimental 
Geometries 
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Figure  3.  Radial  Distribution  Function  (RDF)  for  Mass-Center-Mass-Center  Pairs  as 
Functions  of  Temperature  for  the  e  Phase. 

present,  no  experimental  data  are  available  to  which  the  calculated  thermal  coefficients  can  be 
compared,  and  it  is  hoped  that  these  results  will  stimulate  measurement  of  these  properties. 

3.  Conclusion 

We  have  performed  MP  and  NPT-MD  simulations  of  three  phases  (e,  p,  and  y)  of  the 
2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  (HNIW)  crystal  using  6-exp  Buckingham  potentials 
developed  for  the  RDX  crystal  [1],  plus  Coulombic  interactions  using  electrostatic  charges 
determined  from  fits  to  ab  initio  electrostatic  potentials  calculated  at  the  HF/6-31G**  level.  MP 
calculations  with  and  without  symmetry  constraints  show  good  agreement  between  predicted 
geometrical  parameters  and  experimental  values  for  all  three  phases.  Additionally,  the  calculations 
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Figure  4.  Unit  Cell  Edge  Lengths  and  Volumes  as  Functions  of  Temperature  for  e-HNIW  (Upper  Frames);  p-HNIW  (Middle 
Frames);  and  y-HNIW  (Lower  Frames).  The  Symbols  Denote  Time  Averages,  and  the  Solid  Curves  Are  Linear  Fits  of 
the  Points  (Coefficients  Given  in  Table  6). 
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Table  6.  Calculated  NPT-MD  Lattice  Dimensions  at  Various  Temperatures 


Polymorphic  Phase 

T 

(K) 

a 

(A) 

b 

(A) 

c 

(A) 

Volume 

(A!) 

8.7953 

12.5006 

13.4066 

1422.8883 

100.0 

8.8088 

12.5232 

13.4325 

1429.8001 

200.0 

8.8247 

12.5568 

13.4648 

1439.1609 

273.1 

8.8348 

12.5751 

13.4841 

1444.5670 

e  Phase 

300.0 

8.8420 

12.5837 

13.4897 

1447.1754 

325.0 

8.8446 

12.5903 

13.4959 

1448.8157 

350.0 

8.8467 

12.5976 

13.5048 

1450.8332 

375.0 

8.8548 

12.6056 

13.5128 

1453.9118 

400.0 

8.8552 

12.6083 

13.5182 

1454.5415 

425.0 

8.8598 

12.6189 

13.5298 

1457.2820 

mm 

9.5272 

12.8748 

11.7050 

1435.7761 

100.0 

9.5513 

12.8904 

11.7188 

1442.8274 

200.0 

9.5776 

12.9143 

11.7417 

1452.3071 

273.1 

9.5954 

12.9285 

11.7520 

1457.8772 

P  Phase 

300.0 

9.6105 

12.9316 

11.7566 

1461.0940 

330.0 

9.6116 

12.9385 

11.7653 

1463.11 66 

350.0 

9.6269 

12.9467 

11.7715 

1467.1411 

375.0 

9.6270 

12.9570 

11.7746 

1468.7066 

400.0 

9.6359 

12.9591 

11.7832 

1471.3920 

425.0 

9.6462 

12.9671 

11.7867 

1474.2993 
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Table  6.  Calculated  NPT-MD  Lattice  Dimensions  at  Various  Temperatures  (continued) 


Polymorphic  Phase 

a 

(A) 

b 

(A) 

c 

(A) 

Volume 

(A3) 

y  Phase 

4.2 

13.4370 

7.9092 

14.8548 

1493.9136 

100.0 

13.4535 

7.9281 

14.8786 

1501.3136 

200.0 

13.4771 

7.9490 

14.9036 

1509.8282 

273.1 

13.5038 

7.9643 

14.9318 

1517.8931 

300.0 

13.5144 

7.9690 

14.9413 

1520.8156 

325.0 

13.5205 

7.9740 

14.9510 

1523.1271 

350.0 

13.5232 

7.9755 

14.9510 

1523.9162 

375.0 

13.5296 

7.9840 

14.9635 

1527.2047 

400.0 

13.5515 

7.9868 

14.9705 

1530.8599 

425.0 

13.5500 

7.9915 

14.9824 

1532.3298 

indicate  a  stability  ranking  order  e  >  P  >  y  in  agreement  with  experimental  measurement  [5]. 
Predictions  of  crystal  parameters  at  room  temperature  and  zero  pressure  agree  with  the  experimental 
unit  cell  dimensions  to  within  1.0%  for  phase  e,  0.9%  for  phase  p,  and  2.5%  for  phase  y. 
Additionally,  little  rotational  or  translational  disorder  occurs  in  thermal,  unconstrained  trajectories. 

Temperature  dependencies  of  the  physical  parameters  of  the  lattice  at  zero  pressure  over  the 
temperature  range  4.2-425  K  indicate  that  thermal  expansion  of  the  three  crystalline  polymorphs  is 
nearly  isotropic  for  e-  and  y-HNIW  and  is  anisotropic  for  P-HNIW. 

The  success  of  the  present  potential  energy  parameters  in  describing  different  phases  of  the 
HNIW  crystal  at  moderate  temperatures  and  low  pressures  provides  incentive  to  further  investigate 
the  transferability  of  this  model  to  other  cyclic  nitramine  systems  (e.g.,  HMX)  and  more  dissimilar 
energetic  crystals  (TNT,  TATB).  Future  work  will  be  directed  at  determining  the  degree  to  which 
this  potential  energy  function  is  transferable  to  other  energetic  materials.  Incorporation  of 
intramolecular  motion  by  relaxing  the  rigid  molecular  model  will  also  be  investigated. 
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